Étude de l'erreur due à la diffusion et à la convection dans le cas du SWARM

import sys

sys.path = ['/home/as259691/PycharmProjects/FluidDyn1D'] + sys.path

from src.main import *
from src.plot_fields import *

# %matplotlib inline
rc('figure', figsize=(10,5))
rc('figure', dpi=100)

Paramètres de la simulation

Ici on va réaliser une simulation sans diffusion pour différentes écritures de notre équation thermique.

La résolution se fait à chaque fois en WENO avec Euler explicite en temps.

# d = 6./100*Delta/2.
phy_prop = PhysicalProperties(Delta=0.02, v=0.2, dS=0.005**2,
                              lda1=5.5*10**-2, lda2=15.5, rho_cp1=70278., rho_cp2=702780., diff=1.,
                              alpha=0.06, a_i=357.)
markers = Bulles(phy_prop=phy_prop)
Formulation = [Problem, ProblemConserv2]

n = 1000
t_fin = 0.2

def compare_energy_forme(formu, phy_prop, num_prop, markers, t_fin):
    fig1,ax1 = plt.subplots(1)
    ax1.set_title('Énergie en fonction du temps')
    for form in formu:
        print()
        prob = form(get_T_creneau, markers=markers, phy_prop=phy_prop, num_prop=num_prop)
        E0 = prob.energy
        print(prob.name)
        print('==========================')
        t, e = prob.timestep(t_fin=t_fin, number_of_plots=5, debug=False, plotter=Plotter('decale'))
        l = ax1.plot(t, e/(0.02*0.005*0.005), label=prob.name)
        n = len(e)
        i0 = int(n/5)
        dedt_adim = (e[-1] - e[i0]) / (t[-1] - t[i0]) * prob.dt / E0  # on a mult
        # par Dt / rho_cp_l T_l V
        print('dE*/dt* = %g' % dedt_adim)
    le = fig1.legend()

Étude selon la formulation de l'équation

En fait s'il n'y a pas de convection il n'y a pas de différence entre les différentes formes, à l'exception de la moyenne utilisée pour $\frac{1}{\rho C_p}$

RK4

num_prop = NumericalProperties(dx=3.9*10**-5, schema='weno', time_scheme='rk4', phy_prop=phy_prop, cfl=0.5)
compare_energy_forme(Formulation, phy_prop, num_prop, markers, t_fin)

dt fourier
6.918433404737903e-06
Cas : mixte, rk4, weno, dx = 3.90625e-05, dt = 6.91843e-06, cfl = 0.0354224
==========================
dE*/dt* = -0.000006

dt fourier
6.918433404737903e-06
Forme conservative boniou, Cas : mixte, rk4, weno, dx = 3.90625e-05, dt = 6.91843e-06, cfl = 0.0354224
==========================
dE*/dt* = 0.000022

Euler

num_prop = NumericalProperties(dx=3.9*10**-5, schema='weno', time_scheme='euler', phy_prop=phy_prop)

compare_energy_forme(Formulation, phy_prop, num_prop, markers, t_fin)

dt fourier
6.918433404737903e-06
Cas : mixte, euler, weno, dx = 3.90625e-05, dt = 6.91843e-06, cfl = 0.0354224
==========================
dE*/dt* = -0.000005

dt fourier
6.918433404737903e-06
Forme conservative boniou, Cas : mixte, euler, weno, dx = 3.90625e-05, dt = 6.91843e-06, cfl = 0.0354224
==========================
dE*/dt* = 0.000023

Cas boniou, euler avec weno et upwind quand le stencil traverse l'interface

num_prop = NumericalProperties(dx=3.9*10**-5, schema='weno upwind', time_scheme='euler', phy_prop=phy_prop)
compare_energy_forme(Formulation, phy_prop, num_prop, markers, t_fin)

dt fourier
6.918433404737903e-06
Cas : mixte, euler, weno upwind, dx = 3.90625e-05, dt = 6.91843e-06, cfl = 0.0354224
==========================
dE*/dt* = -0.000006

dt fourier
6.918433404737903e-06
Forme conservative boniou, Cas : mixte, euler, weno upwind, dx = 3.90625e-05, dt = 6.91843e-06, cfl = 0.0354224
==========================
dE*/dt* = 0.000005

Étude numérique

Étude en $\Delta x$

dx_list = [3.9*10**-6, 3.9*10**-5, 1*10**-5]
t_fin = 0.01
for dx in dx_list:
    num_prop = NumericalProperties(dx=dx, schema='weno upwind', time_scheme='euler', phy_prop=phy_prop)
    compare_energy_forme(Formulation, phy_prop, num_prop, markers, t_fin)

dt fourier
6.896863867106643e-08
Cas : mixte, euler, weno upwind, dx = 3.90016e-06, dt = 6.89686e-08, cfl = 0.00353671
==========================
dE*/dt* = -0.000000

dt fourier
6.896863867106643e-08
Forme conservative boniou, Cas : mixte, euler, weno upwind, dx = 3.90016e-06, dt = 6.89686e-08, cfl = 0.00353671
==========================
dE*/dt* = 0.000000

dt fourier
6.918433404737903e-06
Cas : mixte, euler, weno upwind, dx = 3.90625e-05, dt = 6.91843e-06, cfl = 0.0354224
==========================
dE*/dt* = -0.000103

dt fourier
6.918433404737903e-06
Forme conservative boniou, Cas : mixte, euler, weno upwind, dx = 3.90625e-05, dt = 6.91843e-06, cfl = 0.0354224
==========================
dE*/dt* = 0.000034

dt fourier
4.538601983461999e-07
Cas : mixte, euler, weno upwind, dx = 1.0005e-05, dt = 4.5386e-07, cfl = 0.00907267
==========================
dE*/dt* = -0.000004

dt fourier
4.538601983461999e-07
Forme conservative boniou, Cas : mixte, euler, weno upwind, dx = 1.0005e-05, dt = 4.5386e-07, cfl = 0.00907267
==========================
dE*/dt* = 0.000001

dx_list = [7*10**-5]
t_fin = 0.01
for dx in dx_list:
    num_prop = NumericalProperties(dx=dx, schema='weno upwind', time_scheme='euler', phy_prop=phy_prop)
    compare_energy_forme(Formulation, phy_prop, num_prop, markers, t_fin)

dt fourier
2.232841866976439e-05
Cas : mixte, euler, weno upwind, dx = 7.01754e-05, dt = 2.23284e-05, cfl = 0.063636
==========================
dE*/dt* = -0.000469

dt fourier
2.232841866976439e-05
Forme conservative boniou, Cas : mixte, euler, weno upwind, dx = 7.01754e-05, dt = 2.23284e-05, cfl = 0.063636
==========================
dE*/dt* = 0.000172

Étude en $\Delta t$

Ici on ne change pas $\Delta x$, mais on diminue dt_min pour qu'il soit contraignant

dt_list = [2*10**-5, 1*10**-5, 5*10**-6]
t_fin = 0.01
for dt in dt_list:
    num_prop = NumericalProperties(dx=7*10**-5, schema='weno upwind', time_scheme='euler', phy_prop=phy_prop, dt=dt)
    compare_energy_forme(Formulation, phy_prop, num_prop, markers, t_fin)

dt min
2e-05
Cas : mixte, euler, weno upwind, dx = 7.01754e-05, dt = 2e-05, cfl = 0.057
==========================
dE*/dt* = -0.000418

dt min
2e-05
Forme conservative boniou, Cas : mixte, euler, weno upwind, dx = 7.01754e-05, dt = 2e-05, cfl = 0.057
==========================
dE*/dt* = 0.000154

dt min
1e-05
Cas : mixte, euler, weno upwind, dx = 7.01754e-05, dt = 1e-05, cfl = 0.0285
==========================
dE*/dt* = -0.000209

dt min
1e-05
Forme conservative boniou, Cas : mixte, euler, weno upwind, dx = 7.01754e-05, dt = 1e-05, cfl = 0.0285
==========================
dE*/dt* = 0.000076

dt min
4.9999999999999996e-06
Cas : mixte, euler, weno upwind, dx = 7.01754e-05, dt = 5e-06, cfl = 0.01425
==========================
dE*/dt* = -0.000104

dt min
4.9999999999999996e-06
Forme conservative boniou, Cas : mixte, euler, weno upwind, dx = 7.01754e-05, dt = 5e-06, cfl = 0.01425
==========================
dE*/dt* = 0.000038

On n'observe pas de différence notable entre les simulations. Cela laisse à penser que seul le maillage spatial permet de diminuer les pertes d'énergie.

Essayons d'analyser les sources de perte d'énergie.

Étude qualitative

Pertes en diffusion pure

t_fin = 0.01
phy_prop = PhysicalProperties(Delta=0.02, v=0., dS=0.005**2,
                              lda1=5.5*10**-2, lda2=15.5, rho_cp1=70278., rho_cp2=702780., diff=1.,
                              alpha=0.06, a_i=357.)
num_prop = NumericalProperties(dx=7*10**-5, schema='weno upwind', time_scheme='euler', phy_prop=phy_prop)
compare_energy_forme(Formulation, phy_prop, num_prop, markers, t_fin)

dt fourier
2.232841866976439e-05
Cas : diffusion, euler, weno upwind, dx = 7.01754e-05, dt = 2.23284e-05
==========================
dE*/dt* = 0.000010

dt fourier
2.232841866976439e-05
Forme conservative boniou, Cas : diffusion, euler, weno upwind, dx = 7.01754e-05, dt = 2.23284e-05
==========================
dE*/dt* = -0.000000

Perte d'énergie en convection pure

t_fin = 0.01
phy_prop = PhysicalProperties(Delta=0.02, v=0.2, dS=0.005**2,
                              lda1=5.5*10**-2, lda2=15.5, rho_cp1=70278., rho_cp2=702780., diff=0.,
                              alpha=0.06, a_i=357.)
num_prop = NumericalProperties(dx=7*10**-5, schema='weno upwind', time_scheme='euler', phy_prop=phy_prop)
compare_energy_forme(Formulation, phy_prop, num_prop, markers, t_fin)

dt fourier
2.232841866976439e-05
Cas : convection, euler, weno upwind, dx = 7.01754e-05, cfl = 0.063636
==========================
dE*/dt* = -0.000576

dt fourier
2.232841866976439e-05
Forme conservative boniou, Cas : convection, euler, weno upwind, dx = 7.01754e-05, cfl = 0.063636
==========================
dE*/dt* = 0.000166

Étude selon les schémas numériques

Schemas = ['upwind', 'center', 'weno', 'weno upwind']
Time_scheme = ['euler', 'rk4']
def compare_energy_schema(schemas, form, time_scheme, phy_prop, markers, t_fin):
    fig1,ax1 = plt.subplots(1)
    ax1.set_title('Énergie en fonction du temps')
    for schem in schemas:
        for ts in time_scheme:
            a = Plotter('decale')
            num_prop = NumericalProperties(dx=3*10**-5, schema=schem, time_scheme=ts, phy_prop=phy_prop)
            print()
            prob = form(get_T_creneau, markers, phy_prop=phy_prop, num_prop=num_prop)
            E0 = prob.energy
            print(prob.name)
            print('==========================')
            t, e = prob.timestep(t_fin=t_fin, number_of_plots=5, debug=False, plotter=a)
            a.ax.set_xlim(0., phy_prop.Delta/2)
            a.ax.set_ylim(0.6, 1.1)
            l = ax1.plot(t, e/(0.02*0.005*0.005), label=prob.name)
            n = len(e)
            i0 = int(n/5)
            dedt_adim = (e[-1] - e[i0]) / (t[-1] - t[i0]) * prob.dt / E0  # on a mult
            # par Dt / rho_cp_l T_l V
            print('dE*/dt* = %f' % dedt_adim)
    le = fig1.legend()

t_fin = 0.01
fig1, ax1 = plt.subplots(1)
a = Plotter('decale')
num_prop = NumericalProperties(dx=1*10**-5, schema='weno', time_scheme='rk4', phy_prop=phy_prop)
print()
prob_ref = Problem(get_T_creneau, markers, phy_prop=phy_prop, num_prop=num_prop)
E0 = prob_ref.energy
print(prob_ref.name)
print('==========================')
t, e = prob_ref.timestep(t_fin=t_fin, number_of_plots=2, debug=False, plotter=a)
a.ax.set_ylim(0.5,1.05)
a.ax.set_xlim(0., phy_prop.Delta/4.)
l = ax1.plot(t, e/(0.02*0.005*0.005), label=prob_ref.name)
ax1.legend()
ax1.grid(b=True)
n = len(e)
i0 = int(n/5)
dedt_adim = (e[-1] - e[i0]) / (t[-1] - t[i0]) * prob_ref.dt / E0  # on a mult

num_prop = NumericalProperties(dx=5*10**-6, schema='weno', time_scheme='rk4', phy_prop=phy_prop)
print()
prob_ref = Problem(get_T_creneau, markers, phy_prop=phy_prop, num_prop=num_prop)
E0 = prob_ref.energy
print(prob_ref.name)
print('==========================')
t, e = prob_ref.timestep(t_fin=t_fin, number_of_plots=2, debug=False, plotter=a)
a.ax.set_ylim(0.5,1.05)
a.ax.set_xlim(0., phy_prop.Delta/4.)
l = ax1.plot(t, e/(0.02*0.005*0.005), label=prob_ref.name)
ax1.legend()
ax1.grid(b=True)
n = len(e)
i0 = int(n/5)
dedt_adim = (e[-1] - e[i0]) / (t[-1] - t[i0]) * prob_ref.dt / E0  # on a mult
# par Dt / rho_cp_l T_l V
print('dE*/dt* = %f' % dedt_adim)

dt fourier
4.538601983461999e-07
Cas : convection, rk4, weno, dx = 1.0005e-05, cfl = 0.00907267
==========================

dt fourier
1.1335161290322582e-07
Cas : convection, rk4, weno, dx = 5e-06, cfl = 0.00453406
==========================
dE*/dt* = -0.000000

t_fin = 0.01
Schemas = ['upwind', 'weno', 'weno upwind']
compare_energy_schema(Schemas, Problem, Time_scheme, phy_prop, markers, t_fin)

dt fourier
4.0888316389624144e-06
Cas : convection, euler, upwind, dx = 3.003e-05, cfl = 0.0272316
==========================
dE*/dt* = -0.000129

dt fourier
4.0888316389624144e-06
Cas : convection, rk4, upwind, dx = 3.003e-05, cfl = 0.0272316
==========================
dE*/dt* = -0.000130

dt fourier
4.0888316389624144e-06
Cas : convection, euler, weno, dx = 3.003e-05, cfl = 0.0272316
==========================
dE*/dt* = -0.000009

dt fourier
4.0888316389624144e-06
Cas : convection, rk4, weno, dx = 3.003e-05, cfl = 0.0272316
==========================
dE*/dt* = -0.000017

dt fourier
4.0888316389624144e-06
Cas : convection, euler, weno upwind, dx = 3.003e-05, cfl = 0.0272316
==========================
dE*/dt* = -0.000081

dt fourier
4.0888316389624144e-06
Cas : convection, rk4, weno upwind, dx = 3.003e-05, cfl = 0.0272316
==========================
dE*/dt* = -0.000083

Schemas = ['upwind', 'weno', 'weno upwind']
Time_scheme = ['euler', 'rk4']
compare_energy_schema(Schemas, ProblemConserv2, Time_scheme, phy_prop, markers, t_fin)

dt fourier
4.0888316389624144e-06
Forme conservative boniou, Cas : convection, euler, upwind, dx = 3.003e-05, cfl = 0.0272316
==========================
dE*/dt* = 0.000033

dt fourier
4.0888316389624144e-06
Forme conservative boniou, Cas : convection, rk4, upwind, dx = 3.003e-05, cfl = 0.0272316
==========================
dE*/dt* = 0.000032

dt fourier
4.0888316389624144e-06
Forme conservative boniou, Cas : convection, euler, weno, dx = 3.003e-05, cfl = 0.0272316
==========================
dE*/dt* = 0.000011

dt fourier
4.0888316389624144e-06
Forme conservative boniou, Cas : convection, rk4, weno, dx = 3.003e-05, cfl = 0.0272316
==========================
dE*/dt* = 0.000015

dt fourier
4.0888316389624144e-06
Forme conservative boniou, Cas : convection, euler, weno upwind, dx = 3.003e-05, cfl = 0.0272316
==========================
dE*/dt* = 0.000014

dt fourier
4.0888316389624144e-06
Forme conservative boniou, Cas : convection, rk4, weno upwind, dx = 3.003e-05, cfl = 0.0272316
==========================
dE*/dt* = 0.000014

Profiling du code

from line_profiler import LineProfiler
num_prop = NumericalProperties(dx=7*10**-5, schema='weno upwind', time_scheme='rk4', phy_prop=phy_prop, dt=1.)
lp = LineProfiler()
func = [Problem.rk4_timestep, Bulles.indicatrice_liquide, interpolate_form_center_to_face_weno,
       integrale_volume_div]
for fu in func:
    lp.add_function(fu)   # add additional function to profile
lp_wrapper = lp(compare_energy_forme)
lp_wrapper(Formulation, phy_prop, num_prop, markers, 0.1)
lp.print_stats()

dt fourier
2.232841866976439e-05
Cas : convection, rk4, weno upwind, dx = 7.01754e-05, cfl = 0.063636
==========================
dE*/dt* = -0.000018

dt fourier
2.232841866976439e-05
Forme conservative boniou, Cas : convection, rk4, weno upwind, dx = 7.01754e-05, cfl = 0.063636
==========================
dE*/dt* = 0.000018
Timer unit: 1e-06 s

Total time: 13.7396 s
File: /home/as259691/PycharmProjects/FluidDyn1D/src/main.py
Function: integrale_volume_div at line 10

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
    10                                           def integrale_volume_div(center_value, face_value, I=None, cl=1, dS=1., schema='weno'):
    11                                               """
    12                                               Calcule le delta de convection aux bords des cellules
    13                                               :param dS:
    14                                               :param schema:
    15                                               :param center_value: les valeurs présentes aux centres des cellules
    16                                               :param face_value: les valeurs présentes aux faces des cellules
    17                                               :param cl: si cl = 1, on prend des gradients nuls aux bords du domaine
    18                                               :return: le delta au centre
    19                                               """
    20     89560      97023.0      1.1      0.7      if len(center_value.shape) != 1 or len(face_value.shape) != 1:
    21                                                   raise NotImplementedError
    22     89560      64403.0      0.7      0.5      if center_value.shape[0] != face_value.shape[0] - 1:
    23                                                   print('Le tableau des valeurs aux faces ne correspond pas au tableau des valeurs au centre')
    24                                                   print(center_value.shape, face_value.shape)
    25                                                   raise NotImplementedError
    26     89560      35334.0      0.4      0.3      if schema is 'upwind':
    27                                                   interpolated_value = interpolate_from_center_to_face_upwind(center_value, cl)
    28     89560      31055.0      0.3      0.2      elif schema is 'center':
    29                                                   interpolated_value = interpolate_from_center_to_face_center(center_value, cl)
    30     89560      29909.0      0.3      0.2      elif schema is 'weno':
    31                                                   interpolated_value = interpolate_form_center_to_face_weno(center_value, cl=1)
    32     89560      40326.0      0.5      0.3      elif schema == 'weno upwind':
    33     89560      30846.0      0.3      0.2          if I is None:
    34                                                       raise NotImplementedError
    35     89560   13061058.0    145.8     95.1          interpolated_value = interpolate_center_value_weno_to_face_upwind_interface(center_value, I, cl=1)
    36                                               else:
    37                                                   raise NotImplementedError
    38     89560     106096.0      1.2      0.8      flux = interpolated_value * face_value
    39                                               # plt.figure(1)
    40                                               # plt.plot(flux)
    41                                               # plt.show(block=False)
    42     89560     137855.0      1.5      1.0      delta_center_value = flux[1:] - flux[:-1]
    43     89560     105728.0      1.2      0.8      return dS * delta_center_value

Total time: 7.54102 s
File: /home/as259691/PycharmProjects/FluidDyn1D/src/main.py
Function: interpolate_form_center_to_face_weno at line 73

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
    73                                           def interpolate_form_center_to_face_weno(a, cl=1):
    74                                               """
    75                                               Weno scheme
    76                                               :param a: the scalar value at the center of the cell
    77                                               :param cl: conditions aux limites, cl = 1: périodicité
    78                                               :return: les valeurs interpolées aux faces de la face -1/2 à la face n+1/2
    79                                               """
    80     89560     113876.0      1.3      1.5      if cl == 1:
    81     89560     229841.0      2.6      3.0          center_values = np.empty(a.size + 5)
    82     89560     217261.0      2.4      2.9          center_values[:3] = a[-3:]
    83     89560     147485.0      1.6      2.0          center_values[3:-2] = a
    84     89560     153844.0      1.7      2.0          center_values[-2:] = a[:2]
    85                                               else:
    86                                                   raise NotImplementedError
    87     89560     120102.0      1.3      1.6      ujm2 = center_values[:-4]
    88     89560     116337.0      1.3      1.5      ujm1 = center_values[1:-3]
    89     89560     116051.0      1.3      1.5      uj = center_values[2:-2]
    90     89560     118486.0      1.3      1.6      ujp1 = center_values[3:-1]
    91     89560     117561.0      1.3      1.6      ujp2 = center_values[4:]
    92     89560     473192.0      5.3      6.3      f1 = 1. / 3 * ujm2 - 7. / 6 * ujm1 + 11. / 6 * uj
    93     89560     376808.0      4.2      5.0      f2 = -1. / 6 * ujm1 + 5. / 6 * uj + 1. / 3 * ujp1
    94     89560     369830.0      4.1      4.9      f3 = 1. / 3 * uj + 5. / 6 * ujp1 - 1. / 6 * ujp2
    95     89560     150119.0      1.7      2.0      eps = np.array(10. ** -6)
    96     89560     987023.0     11.0     13.1      b1 = 13. / 12 * (ujm2 - 2 * ujm1 + uj) ** 2 + 1. / 4 * (ujm2 - 4 * ujm1 + 3 * uj) ** 2
    97     89560     657626.0      7.3      8.7      b2 = 13. / 12 * (ujm1 - 2 * uj + ujp1) ** 2 + 1. / 4 * (ujm1 - ujp1) ** 2
    98     89560     850713.0      9.5     11.3      b3 = 13. / 12 * (uj - 2 * ujp1 + ujp2) ** 2 + 1. / 4 * (3 * uj - 4 * ujp1 + ujp2) ** 2
    99     89560     334858.0      3.7      4.4      w1 = 1. / 10 / (eps + b1) ** 2
   100     89560     309056.0      3.5      4.1      w2 = 3. / 5 / (eps + b2) ** 2
   101     89560     304737.0      3.4      4.0      w3 = 3. / 10 / (eps + b3) ** 2
   102     89560     207644.0      2.3      2.8      sum_w = w1 + w2 + w3
   103     89560     254552.0      2.8      3.4      w1 /= sum_w
   104     89560     180639.0      2.0      2.4      w2 /= sum_w
   105     89560     171257.0      1.9      2.3      w3 /= sum_w
   106     89560     358700.0      4.0      4.8      interpolated_value = f1 * w1 + f2 * w2 + f3 * w3
   107     89560     103418.0      1.2      1.4      return interpolated_value

Total time: 4.6227 s
File: /home/as259691/PycharmProjects/FluidDyn1D/src/main.py
Function: indicatrice_liquide at line 245

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
   245                                               def indicatrice_liquide(self, x):
   246                                                   """
   247                                                   Calcule l'indicatrice qui correspond au liquide avec les marqueurs selon la grille x
   248                                                   :param x: les positions des centres des mailles
   249                                                   :return: l'indicatrice
   250                                                   """
   251     44784     418997.0      9.4      9.1          i = np.ones_like(x)
   252     44784      56580.0      1.3      1.2          dx = x[1] - x[0]
   253    223920     220913.0      1.0      4.8          for markers in self.markers:
   254    179136     131025.0      0.7      2.8              if markers[0] < markers[1]:
   255    176447     778653.0      4.4     16.8                  i[(x > markers[0]) & (x < markers[1])] = 0.
   256                                                       else:
   257      2689      12116.0      4.5      0.3                  i[(x > markers[0]) | (x < markers[1])] = 0.
   258    179136     687659.0      3.8     14.9              diph0 = (np.abs(x - markers[0]) < dx / 2.)
   259    179136     834643.0      4.7     18.1              i[diph0] = (markers[0] - x[diph0]) / dx + 1. / 2.
   260    179136     630356.0      3.5     13.6              diph1 = (np.abs(x - markers[1]) < dx / 2.)
   261    179136     833671.0      4.7     18.0              i[diph1] = -(markers[1] - x[diph1]) / dx + 1 / 2.
   262     44784      18085.0      0.4      0.4          return i

Total time: 10.1102 s
File: /home/as259691/PycharmProjects/FluidDyn1D/src/main.py
Function: rk4_timestep at line 528

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
   528                                               def rk4_timestep(self, debug=False):
   529      4478      11374.0      2.5      0.1          T_int = self.T.copy()
   530      4478       5011.0      1.1      0.0          K = [0.]
   531      4478      17374.0      3.9      0.2          pas_de_temps = np.array([0, 0.5, 0.5, 1.])
   532     22390      34396.0      1.5      0.3          for h in pas_de_temps:
   533     17912     346440.0     19.3      3.4              markers_int = self.markers.copy()
   534     17912     179447.0     10.0      1.8              markers_int.shift(self.phy_prop.v * h * self.dt)
   535     17912    2151399.0    120.1     21.3              temp_I = markers_int.indicatrice_liquide(self.num_prop.x)
   536     17912      72496.0      4.0      0.7              T = T_int + h * self.dt * K[-1]
   537     17912      36542.0      2.0      0.4              int_div_T_u = -1 / (self.phy_prop.dS * self.num_prop.dx) * \
   538     17912     145813.0      8.1      1.4                  integrale_volume_div(T, self.phy_prop.v * np.ones((T.shape[0] + 1,)), I=temp_I,
   539     17912    2977345.0    166.2     29.4                                       dS=self.phy_prop.dS, schema=self.num_prop.schema)
   540     17912     127346.0      7.1      1.3              Lda_h = 1. / (temp_I / self.phy_prop.lda1 + (1. - temp_I) / self.phy_prop.lda2)
   541     17912      86600.0      4.8      0.9              rho_cp_inv_h = temp_I / self.phy_prop.rho_cp1 + (1. - temp_I) / self.phy_prop.rho_cp2
   542     17912      31650.0      1.8      0.3              div_lda_grad_T = 1 / (self.phy_prop.dS * self.num_prop.dx) * \
   543     17912     187282.0     10.5      1.9                  integrale_volume_div(Lda_h, grad(T, dx=self.num_prop.dx), I=temp_I,
   544     17912    2860609.0    159.7     28.3                                       dS=self.phy_prop.dS, schema=self.num_prop.schema)
   545     17912      58233.0      3.3      0.6              rho_cp_inv_int_div_lda_grad_T = self.phy_prop.diff * rho_cp_inv_h * div_lda_grad_T
   546     17912      33978.0      1.9      0.3              K.append(int_div_T_u + rho_cp_inv_int_div_lda_grad_T)
   547     17912      16489.0      0.9      0.2              if debug:
   548                                                           plt.figure('sous-pas de temps %f' % (len(K) - 2))
   549                                                           plt.plot(self.num_prop.x_f, interpolate_form_center_to_face_weno(Lda_h) * grad(T, dx=self.num_prop.dx),
   550                                                                    label='lda_h grad T, time = %f' % self.time)
   551                                                           plt.plot(self.num_prop.x, rho_cp_inv_h, label='rho_cp_inv_h, time = %f' % self.time)
   552                                                           plt.plot(self.num_prop.x, div_lda_grad_T, label='div_lda_grad_T, time = %f' % self.time)
   553                                                           maxi = max(np.max(div_lda_grad_T), np.max(rho_cp_inv_int_div_lda_grad_T), np.max(rho_cp_inv_h))
   554                                                           mini = min(np.min(div_lda_grad_T), np.min(rho_cp_inv_int_div_lda_grad_T), np.min(rho_cp_inv_h))
   555                                                           for markers in self.markers():
   556                                                               plt.plot([markers[0] + self.phy_prop.v * h * self.dt] * 2, [mini, maxi], '--')
   557                                                               plt.plot([markers[1] + self.phy_prop.v * h * self.dt] * 2, [mini, maxi], '--')
   558                                                               plt.xticks(self.num_prop.x_f)
   559                                                               plt.grid(b=True, which='major')
   560                                                               plt.legend()
   561      4478      17144.0      3.8      0.2          coeff = np.array([1. / 6, 1 / 3., 1 / 3., 1. / 6])
   562      4478     116848.0     26.1      1.2          self.T += np.sum(self.dt * coeff * np.array(K[1:]).T, axis=-1)
   563      4478     596414.0    133.2      5.9          self.update_markers()

Total time: 27.6902 s
File: <ipython-input-4-6361008d467a>
Function: compare_energy_forme at line 1

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================

Comparaison des meilleurs cas pour chaque formulation

t_fin = 0.2
fig1, ax1 = plt.subplots(1)

for dx in [2*10**-5, 4*10**-5, 7*10**-5, 10*10**-5]:
    print('~~~~~~~~~~~~~~~~~~~~~~~~~~')
    a = Plotter('decale')
    num_prop = NumericalProperties(dx=dx, schema='weno', time_scheme='rk4', phy_prop=phy_prop)
    print()
    prob_ref = Problem(get_T_creneau, markers, phy_prop=phy_prop, num_prop=num_prop)
    E0 = prob_ref.energy_m
    print(prob_ref.name)
    print('==========================')
    t, e = prob_ref.timestep(t_fin=t_fin, number_of_plots=2, debug=False, plotter=a)
    l = ax1.plot(t, e, label=prob_ref.name)
    n = len(e)
    i0 = int(n/5)
    dedt_adim = (e[-1] - e[i0]) / (t[-1] - t[i0]) * prob_ref.dt / E0  # on a mult
    print('dE/dt = %g' % dedt_adim)

    num_prop = NumericalProperties(dx=dx, schema='weno upwind', time_scheme='rk4', phy_prop=phy_prop)
    print()
    prob_ref = ProblemConserv2(get_T_creneau, markers, phy_prop=phy_prop, num_prop=num_prop)
    E0 = prob_ref.energy_m
    print(prob_ref.name)
    print('==========================')
    t, e = prob_ref.timestep(t_fin=t_fin, number_of_plots=2, debug=False, plotter=a)
    l = ax1.plot(t, e, label=prob_ref.name)
    n = len(e)
    i0 = int(n/5)
    dedt_adim = (e[-1] - e[i0]) / (t[-1] - t[i0]) * prob_ref.dt 
    print('dE/dt = %g' % dedt_adim)
    
    a.ax.set_ylim(0.,1.05)
    a.ax.set_xlim(0., phy_prop.Delta/4.)
    ax1.legend()
    ax1.grid(b=True)

~~~~~~~~~~~~~~~~~~~~~~~~~~

dt fourier
1.8172585062055175e-06
Cas : convection, rk4, weno, dx = 2.002e-05, cfl = 0.0181544
==========================
dE/dt = -1.69281e-13

dt fourier
1.8172585062055175e-06
Forme conservative boniou, Cas : convection, rk4, weno upwind, dx = 2.002e-05, cfl = 0.0181544
==========================
dE/dt = 5.65927e-09
~~~~~~~~~~~~~~~~~~~~~~~~~~

dt fourier
7.283608525474247e-06
Cas : convection, rk4, weno, dx = 4.00802e-05, cfl = 0.0363452
==========================
dE/dt = -2.01575e-12

dt fourier
7.283608525474247e-06
Forme conservative boniou, Cas : convection, rk4, weno upwind, dx = 4.00802e-05, cfl = 0.0363452
==========================
dE/dt = 3.91957e-08
~~~~~~~~~~~~~~~~~~~~~~~~~~

dt fourier
2.232841866976439e-05
Cas : convection, rk4, weno, dx = 7.01754e-05, cfl = 0.063636
==========================
dE/dt = -8.5763e-12

dt fourier
2.232841866976439e-05
Forme conservative boniou, Cas : convection, rk4, weno upwind, dx = 7.01754e-05, cfl = 0.063636
==========================
dE/dt = 1.82845e-07
~~~~~~~~~~~~~~~~~~~~~~~~~~

dt fourier
4.534064516129032e-05
Cas : convection, rk4, weno, dx = 0.0001, cfl = 0.0906813
==========================
dE/dt = -1.45701e-11

dt fourier
4.534064516129032e-05
Forme conservative boniou, Cas : convection, rk4, weno upwind, dx = 0.0001, cfl = 0.0906813
==========================
dE/dt = 4.96835e-07

markers.shift(10**-6/7.)

%matplotlib inline
rc('figure', figsize=(10,7))
rc('figure', dpi=100)
n = 1000
t_fin = 10**-5
fig1, ax1 = plt.subplots(1)


for dx in [5*10**-6, 1*10**-6, 5*10**-7]:
    print('~~~~~~~~~~~~~~~~~~~~~~~~~~')
    a = Plotter('decale')
    num_prop = NumericalProperties(dx=dx, schema='weno', time_scheme='rk4', phy_prop=phy_prop)
    print()
    prob_ref = Problem(get_T_creneau, markers=markers, phy_prop=phy_prop, num_prop=num_prop)
    E0 = prob_ref.energy
    print(prob_ref.name)
    print('==========================')
    t, e = prob_ref.timestep(n=n, t_fin=t_fin, number_of_plots=1, plotter=a)
    l = ax1.plot(t, e/(0.02*0.005*0.005), label=prob_ref.name)
    n = len(e)
    i0 = int(n/5)
    dedt_adim = (e[-1] - e[i0]) / (t[-1] - t[i0]) * prob_ref.dt / E0  # on a mult
    print('dE/dt = %g' % dedt_adim)

    num_prop = NumericalProperties(dx=dx, schema='weno upwind', time_scheme='rk4', phy_prop=phy_prop)
    print()
    prob_ref = ProblemConserv2(get_T_creneau, markers=markers, phy_prop=phy_prop, num_prop=num_prop)
    E0 = prob_ref.energy
    print(prob_ref.name)
    print('==========================')
    t, e = prob_ref.timestep(t_fin=t_fin, n=n, number_of_plots=1, plotter=a)
    l = ax1.plot(t, e/(0.02*0.005*0.005), label=prob_ref.name)
    n = len(e)
    i0 = int(n/5)
    dedt_adim = (e[-1] - e[i0]) / (t[-1] - t[i0]) * prob_ref.dt 
    print('dE/dt = %g' % dedt_adim)
    
    a.ax.set_ylim(0.9,1.05)
    a.ax.set_xlim(prob_ref.markers.markers[0][0] - prob_ref.markers.diam / 4., prob_ref.markers.markers[0][0] + prob_ref.markers.diam / 2.,)
    ax1.legend()
    ax1.grid(b=True, which='both')

~~~~~~~~~~~~~~~~~~~~~~~~~~

dt fourier
1.1335161290322582e-07
Cas : mixte, rk4, weno, dx = 5e-06, dt = 1.13352e-07, cfl = 0.00453406
==========================
dE/dt = -2.04651e-05

dt fourier
1.1335161290322582e-07
Forme conservative boniou, Cas : mixte, rk4, weno upwind, dx = 5e-06, dt = 1.13352e-07, cfl = 0.00453406
==========================
dE/dt = 2.13251e-07
~~~~~~~~~~~~~~~~~~~~~~~~~~

dt fourier
4.534064516129032e-09
Cas : mixte, rk4, weno, dx = 1e-06, dt = 4.53406e-09, cfl = 0.000906813
==========================
dE/dt = -3.07688e-08

dt fourier
4.534064516129032e-09
Forme conservative boniou, Cas : mixte, rk4, weno upwind, dx = 1e-06, dt = 4.53406e-09, cfl = 0.000906813
==========================
dE/dt = 6.01558e-09
~~~~~~~~~~~~~~~~~~~~~~~~~~

dt fourier
1.133516129032258e-09
Cas : mixte, rk4, weno, dx = 5e-07, dt = 1.13352e-09, cfl = 0.000453406
==========================
dE/dt = -8.76081e-09

dt fourier
1.133516129032258e-09
Forme conservative boniou, Cas : mixte, rk4, weno upwind, dx = 5e-07, dt = 1.13352e-09, cfl = 0.000453406
==========================
dE/dt = 7.62287e-10